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Redshift-space distortions in galaxy snrveys happen along the radial direction, breaking statistical 
translation invariance. We constrnct estimators for radial distortions that, using only Fast Fourier 
Transforms (FFTs) of the overdensity field multipoles for a given survey geometry, compute the 
power spectrum monopole, quadrupole and hexadecapole, and generalize snch estimators to the 
bispectrnm. Using realistic mock catalogs we compare the signal to noise of two estimators for the 
power spectrum hexadecapole that require different nnmber of FFTs and measnre the bispectrum 
monopole, qnadrupole and hexadecapole. The resnlting algorithm is very efficient, e.g. for the BOSS 
survey requires about three minutes for ^ = 0, 2, 4 power spectra for scales up to A: = 0.3/iMpc“^ 
and about fifteen additional minutes for £ = 0,2,4 bispectra for all scales and triangle shapes up 
to fe = 0.2 fiMpc”^ on a single core. The speed of these estimators is essential as it makes possible 
to compute covariance matrices from large number of realizations of mock catalogs with realistic 
survey characteristics, and paves the way for improved constrains of gravity on cosmological scales, 
inflation and galaxy bias. 


I. INTRODUCTION 


Redshift-space distortions mm of galaxy clustering are 
key in understanding the three-dimensional distribution 
of large-scale structure and are also a major probe for 
constraining gravity on cosmological scales, as evidenced 
in recent work [3Hin]. Such distortions change the Fourier 
modes from their undistorted real-space values depend¬ 
ing on the orientation of the wave-vectors with respect to 
the line of sight. In modern era surveys with large solid 
angles the line of sight is significantly space-dependent 
(as the radial direction varies over the sky), which makes 
Fourier-space analysis non-trivial beyond the lowest mul¬ 
tipole, the monopole. 

To see this we recall that redshift-space positions s are 
given in terms of real-space positions x by 

s = X — f X {u- x) (1) 

where / = d In jd In a in terms of the linear growth 
factor ZJ+ and scale factor a, and the peculiar velocity 
V = —Tifu with PL = din a/dr the comoving Hubble 
constant (and r conformal time). This means that in 
linear perturbation theory the redshift-space density fluc¬ 
tuations are given by 


(5s(x) = d(x) -b / V • 


x{vl- x) 


( 2 ) 


or ds = d + / (2 Ur/r + drUr) which when the solid angle 
of the survey is small enough (the so-called plane-parallel 
approximation, x —)■ z), goes to ds(x) = d(x) + 
leading to ds(k) = (1 -|- //i^)d(k) in Fourier space [T]. 
Using that in linear theory d = V • u, we can also write 


ds(x) = 1 + f ^ d(x) = d(x) (3) 
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which gives the form of the linear redshift-space distor¬ 
tion operator Vg. A fundamental property of radial dis¬ 
tortions is that, unlike the plane-parallel case, the reshift- 
space map does not commute with translations (as the 
observer defines a privileged location) and thus Vg is not 
an eigenfunction of plane-waves, the effect of distortions 
on a mode is not an eigenvalue (1 -b //i^) anymore, and 
as a result the Fourier power spectrum is no longer diag¬ 
onal m- That is, 

(d«(ki)d,(k2)) =P(ki,k2) (4) 

with P(ki,k 2 ) only becoming diagonal in the plane- 
parallel limit, when P(ki,k 2 ) —t P(ki)dD(ki 2 ) with 
ki 2 = ki -b k 2 . The fact that P(ki,k 2 ) is not diagonal 
is directly related to the space-dependent unit vectors in 
Eq. (§ which makes Dg in Fourier space an integral op¬ 
erator inducing mode-coupling even in linear theory [12j . 
While in principle such matrix contains all the informa¬ 
tion it is not clear how to make simple use of it. One 
would like to condense all the information into a set of 
multipoles as a function of a scalar k as in the plane- 
parallel limit but there appears no simple way to do so. 

Of course, radial distortions preserve isotropy about 
the observer and thus spherical harmonics become a nat¬ 
ural basis for angular modes. Spherical harmonics trans¬ 
form alternatives to Fourier analysis exist and are well 
known, at least for the power spectrum (e.g. [n na¬ 
ns]), although their application to surveys is not done as 
often due to their computational cost. For this reason in 
this paper we concentrate on Fourier analysis and how to 
tackle fast estimation of redshift-space power spectrum 
and bispectrum multipoles in the presence of radial dis¬ 
tortions. 

The plan for the paper is as follows. In section |II A| 
we now discuss a slight generalization of the power spec¬ 
trum and bispectrum estimators for “local” regions in 
space where the fluctuations can be taken to be approx¬ 
imately statistically homogeneous (and thus these local 
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statistics can be taken as diagonal). From these in sec¬ 
tions IIB and IIC we build multipoles estimators for the 
power spectrum and bispectrum respectively, as if each 
of these regions is in the plane-parallel approximation, 
giving us estimators that apply to the general case of ra¬ 
dial distortions. In section |Tn] we discuss the application 
to galaxy surveys and in section IIV] we conclude. 


II. DISTORTIONS ESTIMATORS 
A. Local Estimators 

Since in the presence of radial redshift-space distor¬ 
tions the resulting redshift-space density field is no longer 
statistically homogeneous, it makes sense to define a local 
spectrum density estimator at x, 

.Piocai(k,x) = j Ss{yi + wi)Ss{x + vf 2 ) 

( 5 ) 

where the are the coordinates of the x^ from the center 
of mass of the pair x = (xi -|- X2)/2, that is x^ = x -|- 
and xi2 = xi — X2. Therefore wi = —X12/2, W2 = 
X12/2. Equation © is the Fourier transform of the local 
contribution at x to the correlation function at separation 
X12. The local power spectrum density is real but not 
positive definite. 

We can write Eq. (§ in terms of Fourier coefficients, 

Piocai(k,x) = J 5 s(k-bq/ 2 )( 5 s(-k-bq/ 2 ) (6) 

Integrating over space Eq. (H over space the local power 
density is simply 

y ^^^-Piocai(k,x) = |d,(k )|2 ( 7 ) 

that is, the standard power spectrum estimator when it 
makes sense to average over space (when translation in¬ 
variance holds). Equation ([^ has an expectation value, 

(Piocai(k, x)) = J d\ P(k + q/2, -k + q/2) e**! " (8) 

where we have used that the power spectrum in the pres¬ 
ence of radial distortions is no longer diagonal due to the 
loss of statistical homogeneity or translation invariance. 
In the case that there is statistical translation invariance 
(e.g. in the plane-parallel limit), P(k-|-q/2, —k-|-q/2) = 
P(k)( 5 D(q) and we have that (Piocai(k,x)) = P(k), as 
expected. But note that in general, (Piocai(k, x)) con¬ 
tains the same information as the matrix P(ki,k2), as 
the latter can be recovered from the former by a Fourier 
transform. The advantage of the former now becomes 
obvious, as there is a trace of real space (and thus line of 
sight) that can be used to take multipoles and thus com¬ 
press the information on redshift distortions in a similar 


way as in the plane-parallel limit, as we discuss in the 
next section. 

We can now extend these definitions to write down a 
local bispectrum estimator (ks = —ki2) 




qki.xi 3 +k..x, 3 ) 

(27r)3 (27r)3 


3 

X i5s(x-b Wi) 
2=1 


( 9 ) 


where the Wi are the vectors from the centroid of the 
triangle x to the vertices Xi, i.e. x = (xi - 1 - X2 -I- X 3)/3 
with Xi = X -|- Wi and x^ = Xi — Xj. The centroid 
coordinates of the vertices are wi = 2x13/8 — X23/3, 
W2 = 2x23/3 - X13/3, and W3 = -X13/3 - X23/3. In 
terms of Fourier coefficients. 


5 l° 2 f(x) 




q e 


— zqx 


3 

^^(kj -bq/3) 


( 10 ) 


with expectation value 


(Pl°2/r'(x)) = y d^gP(ki-hq/3,k2-bq/3,k3-bq/3)e 

( 11 ) 

which involves the (non-diagonal) bispectrum (note the 
similarity with Eq. . Again, as in the power spectrum 
case, when it does make sense to spatially average, we 
have 

= ( 12 ) 

the standard bispectrum estimator for a statistically ho¬ 
mogeneous field. 


B. Power Spectrum Multipoles 


From Eq. ([^ we can build the natural estimator of 
power multipoles by using x as the line of sight direction 
to the pair of points. 


Pe{k) = { 2 i+ 1 ) J ^ J Piocai(k,x)/:^(fc • i) 

( 13 ) 

where by angular integration we actually mean integra¬ 
tion over a thin shell in Eourier space centered at fc, for 
any function F 



(14) 


where Nk = Jf. d^q = Airk'^Sk is the volume of the shell 
in fc-space, with Sk the bin size. Equation ( 13 ) can be 
rewritten using Eq. ([^ as, 


Piik) = { 2 £ + l) 


47r J (27r)3 (27r)3 


X (5s(xi)5s(x2) 


(15) 
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where y = xi — X 2 and x = Xi + X 2 . This natural es¬ 
timator is for £ = 2 precisely the same as the so-called 
Yamamoto estimator [16j . The complication with such 
estimators is well-known, i.e. they are expensive to com¬ 
pute beyond £ = 0 because the integrals in Eq. (15) do 


not decouple into a product of Fourier transforms due to 
the X dependence in Ci{k ■ x), and thus when discretized 
the integrals become double sums that are quadratic in 
the number of grid points (or galaxies), leading to an 
bottleneck. 

Given this complication, it has been proposed [Bin] 
to make the replacement for the line of sight definition. 


jC{(k ■ x) Ci{k ■ xi) 


(16) 


to make the two integrals (or sums in the discrete case) 
factorize. Still, without further treatment, one the in¬ 
tegrals (or sums) which contains the Legendre polyno¬ 
mial is not of Fourier form due to the k dependence for 
£ > 0. While faster than the natural estimator which re¬ 
quires dealing with pairs of points, this is still expensive 
to compute compared to the power spectrum monopole 
which can just be computed with a single Fourier trans¬ 
form. The replacement in Eq. (16) has recently been 


found to be rather accurate compared to the natural es¬ 
timator by |18j . and the latter accurate compared to the 
full description of the power spectrum including wide- 
angle effects m, so it is of importance to find a fast way 
to compute it. 

One of the main points of this paper is to point out 
that while the replacement in Eq. (16) cannot be writ¬ 


ten as a product of Fourier transforms, it can in fact be 
written as the sum of a product of Fourier transforms, 
and as a result of this, estimators can be built that are 
computable using a handful of FFTs. Indeed, as a result 
of this replacement, Eq. (15) becomes the cross corre¬ 


lation of the local multipole overdensity with the local 
monopole [TumiiH] where 


<5^(k) = 


(Px 

(^ 


....—ik-: 


Ss{x.)Ci{k ■ x) (17) 


which in fact can be computed by FFTs by simply fac¬ 
torizing out the fc-dependence, e.g. for £ = 2 


with 


<52(k) = ^kikjQ,j(k) - ^5o(k) 


fti(k) = J **"'^5s(x)ii% 


(18) 


(19) 


which depends on volume shape through the cosines Xi 
(as the integral in Eq. 19 for a survey becomes over the 
region where 6s is observed, see section [Hi] below). Since 
this is a symmetric tensor only 6 FFTs are needed to 
compute it in the absence of any symmetry of the survey 


geometry, i.e. 

3 r 


^2 = '^\klQ. 


^"yQyy + k^Qz 
~\~^kxkyQxy H” ‘^kykzQyz ‘^kz^xQz 




( 20 ) 


In the plane-parallel limit, all Qij vanish but Qzz and we 
recover the standard plane-parallel results. For £ = 4 we 
have, similarly 

35 * - - - 5 7 

54(k) = —hkjkikkQ^jikik) - -52(k) - -5o(k) (21) 

o z o 


with 




QijlkO^) — J ^27r)3 ^ 6s{pc) XiXjXiX}^ 


( 22 ) 


and so on. Since Qijik is a fully symmetric tensor, in 
the absence of any symmetries, one needs to compute 15 
FFTs to fully characterize it. 


(54 = 


35 


4' 

.2 £.2, 


Q xxxx 


+ (3) eye. -b 4k-^ky Q^^xxy + (6) eye. 
+Qklky Qxxyy + (3) eye. -b Uklkykz Q xxyz 
-b(3) eye. - ^2 - ^<5o 

(23) 


where the number in parenthesis denotes how many total 
terms belong to each cyclic permutation (and counts the 
numbers of FFTs needed). Again, in the plane-parallel 
limit, all Qiju vanish but Qzzzz and we recover the stan¬ 
dard plane-parallel results. 

The power spectrum multipoles estimator then would 
be 


Pi{k) = (2£ + 1) y ^ 5,(k) 5o(-k) (24) 

and thus for £ = 0, 2,4 a total of 1 -b 6 -b 15 = 22 FFTs 
would be needed. Because of A^logiV scaling even this 
many FFTs is still many orders of magnitude faster than 
the naive N'^ procedure where the Legendre polynomials 
are not written in factorized form. 

One obvious question that arises is whether one can do 
better for the hexadecapole (and higher multipoles) as it 
is rather disappointing that one needs so much additional 
computational cost (an additional 15 FFTs for £ = 4) to 
describe a quantity that has rather poor signal to noise 
compared to £ = 0, 2 at large scales. The answer is that 
one can do better by not constraining oneself to a single 
line of sight as £ > 2 is considered. To see this let us 
consider the case £ = 4 for definiteness. We go back to the 
natural estimator in Eq. (15) and write the fourth-order 


Legendre polynomial in terms of quadratic combinations 
of lower even multipoles. 
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U{k-x) = ^ 


C 2 {k ■ x) 


(25) 


where the indices 1,2,3 in ^ow refer to the q^, which 
are being averaged over a shell of thickness 5k about the 
ki, i.e. (pQi = J for thin shells, and 


and then we simply split the first term 


C 2 {k-x) C 2 {k ■ Xi) C 2 {k ■ X 2 ) (26) 


leading to the factorized hexadecapole estimator, 

A 6 (fc) = y f ^|(52(k)|2-P2(fc)-Ipo(fc) (27) 


which does not require any additional FFTs over those 
already computed for £ = 0 , 2 and is based on the auto¬ 
correlation of the local quadrupole of overdensities gen¬ 
erated by the redshift-space mapping. As a result of this 
split, computing £ = 0,2,4 requires only 7 FFTs instead 
of 22 , leading to additional computational savings of just 
over a factor of three by using P 41 , instead of P 4 . In [18] 
it was found that the estimator P 4 in Eq. (241 has a 


small bias compared to the natural estimator at large 
scales, we study the bias o f Pt h relative to P 4 and their 
cosmic variance in section 


III 


below, and find that P 4 
is preferred due to lower cosmic variance, although the 
difference might be negligible for future surveys. 

The same line-of-sight split trick can be used for higher 
multipoles, in the obvious way. Note that for £ = 6 one 
cannot avoid computing the 15 FFTs, but this also allows 
one to compute £ = 8 without extra FFTs. In general 
each new set of FFTs becomes useful for two multipoles, 
as the Legendre polynomials can be split in quadratic 
combinations because of the two line of sights available 
in a two-point function. 


C. Bispectrum Multipoles 

We now consider the case of the bispectrum. In the 
plane-parallel limit, the bispectrum becomes a function of 
five variables: the three sides plus two angular variables 
describing the orientation of the triangle with respect to 
the line of sight. A third angular variable is irrelevant 
in the sense that it rotates the triangle about the line of 
sight, leaving the redshift-space bispectrum invariant. A 
convenient way to handle the plane-parallel case is then 
to do a spherical harmonic decomposition with respect 
to the two relevant angular variables. Let ki > k 2 > k^ 
without loss of generality. From the local bispectrum 
estimator in Eq. (|^ we can then define the multipoles in 
the radial distortions case by 

Bi23 = —11/ '^qifc(qi23) 

’ 123 i—\ 

^Bl°f(x) W„(01,</>12) (28) 


A^123 = 


^ r 

/ (i^gifc(qi23) - 87r^fci^2fc3<5fc^ (29) 

j_l Jki 


In Eq. (28 1 we followed [20] where cos 9i = qi-x, cos 642 = 


qi ■ q 2 and (f>i 2 is the azimuthal angle of q 2 around qi sat¬ 
isfying cos 62 = cos 01 cos 012 — sin0i sin 0 i 2 cos (j)i 2 . An¬ 
other possible choice of angular variables are those that 
describe the orientation of the normal to the triangle face. 
However, it has the disadvantage that it is not well de¬ 
fined for zero area triangles but this can be handled sep¬ 
arately as for such triangles a Legendre multipole decom¬ 
position is all that is needed as all qi differ at most by a 
sign (irrelevant for even multipoles). 

Here for simplicity we take Legendre multipoles with 
respect to the largest side, corresponding to to = 0 mul¬ 
tipoles in Eq. (28 1 , that is 




-^1^3 




J (27r)3 

which can also be written as. 


B1°THx) Ciiqi ■ x) (30) 


-''123 dki 

X jkliiqi ■ x) (31) 

which is obviously in the same “natural” form as for the 
power spectrum. We must now deal with the separability 
of the estimator, as in the power spectrum case. For £ = 2 
using X ^ xi we obtain the factorized estimator for the 
bispectrum quadrupole. 


^123 = 5 n / ^P(^23) ^^(qi) ,5o(q2) ,5o(q3) 

i=l -''123 

(32) 

and so computing the bispectrum quadrupole only results 
in a total factor of two over monopole alone, as the com¬ 
putation of 62 is negligible in cost with the bispectrum 
itself (and is the same ingredient needed for the power 
spectrum quadrupole). A similar consideration leads to 
the hexadecapole bispectrum estimator. 


-§123 = 9 n / ^^(qg) 

i=l -''123 

(33) 

and the obvious generalization for higher-order multi¬ 
poles. As discussed above for the power spectrum, addi¬ 
tional multipole information may be obtained instead by 
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including more than one quadrupole field 82 in Eq. (32) 
instead of computing the additional 15 FFTs to build 84 ^. 
That is, if for £ = 4 we may use Eq. (31) and Eq. (25) 
with the split 


C2{qi-x) C2{qi ■ xi) C2{qi ■ X 2 ) 


(34) 


we obtain the alternative hexadecapole bispectrum esti¬ 
mator, 


B 


(4b) 

123 


tII / ^2(qi)^2(9i,q2)i5o(q3) 

^ i = l -'''123 

(35) 


where 


3 1 

^ 2 (p, q) = 2 ‘ 3 b (q) - 2'^o('^) (^®) 

Note that for zero area triangles < 52 ( 91 , q 2 ) = <52(q2) and 
thus, 


rNI = 

^123 — 


35 


n 


[ rf^q* ^2(qi) ^2(q2) ^o(q3) 

Iki -‘''n- 

7 


'123 


_5(2)_^(0) 

-°123 2 “°123 


(37) 


which is analogous to the power spectrum case Eq. (27). 


A disadvantage of the estimator in Eq. (35) is that the 
extra qi dependence means that the bispectrum estima¬ 
tor is a bit more costly, as e.g. for £ = 4 we must now 
estimate 


/ d^qifc(qi23) (9i)*(9i)j i52(qi)Q*j(q2)i5o(q3) 

i=l 

(38) 

which corresponds to evaluating six bispectra. Since the 
cost of evaluating the bispectrum is much more than eval¬ 
uating the next set of 15 FFTs of the overdensity fields 
in Eq. (23), it is more convenient in this case to avoid the 


split in Eq. (34) and instead use Eq. (33). 


Another approach would be to use the zero-area trian¬ 


gle estimator in Eq. (37) for all triangles. Unlike Eq. (33), 


this does not need the extra 15 FFTs and it is as fast to 
estimate. However, this estimator for general triangles 
is not a true multipole, that is, it does not vanish in 
real space except for zero area triangles and therefore ex¬ 
tracting information from it about redshift-space distor¬ 
tions may be more complicated due to degeneracies with 
the monopole. However, it may be useful not just for 
zero area triangles, but at bit more generally for nearly 
squeezed or folded triangles. It is beyond the scope of 
this paper to explore this further, in what follows we will 
consider the more standard bispectrum multipoles as in 
Eq. (32) and (33), particularly in light of the reduced cos¬ 


mic variance of such estimators compared to those that 


require less FFTs as we find below for the power spec¬ 
trum multipoles. 

Finally, for completeness we briefly mention how to 
build the zero-area triangle estimator for £ = 6. We split 
the sixth-order Legendre polynomial in cubic combina¬ 
tions of lower-order even polynomials, as we now have 
three lines of sight. Using that 


Ceik-x) = ^ C2{k-x) 
io L 


we simply split the first term 


-13 7 r ^ 127 . 2 

- C 2 {k-x) --C 2 {k-x) + - 

3 L Jo 9 

(39) 


3 


^ 2(91 • x) 


£2(91 • Xi) £2(92 • X2) £2(93 • <^3) 


(40) 


which leads to, 


= 

^123 — 


1001 


n/ ^^qi '^2(qi)<52(q2)l52(q3) 

Jki -'''123 


26 

Y^^123 


““ r(2) ^0-3 S(0) 


10' 


'123 


283 


45 


123 


(41) 


so now with the six FFTs needed for the quadrupole we 
can compute up to the £ = 6 zero-area triangle bispec¬ 
trum multipole, so each set of FFTs is enough for three 
multipoles as there are three lines of sight in a three-point 
function. 


III. IMPLEMENTATION IN GALAXY 
SURVEYS 


A. Power Spectrum 


We now proceed to implementing the above ideas in 
the case of a survey geometry, where the galaxy sample 
is characterized by Ng galaxies at positions and the 
radial and angular selection functions by a random cat¬ 
alog with Nr objects {a = Ng/Nr 1). Each object is 
given a weight wj, e.g. the FKP weights [2T]. Given the 
results above, we write the overdensity monopole 

Ng Nr 

Fo(k)^(^-a^)u;,e*‘^-' (42) 

1=1 1=1 

whereas for the quadrupole we have 

F2(k) = hahQ^^ik) - ^Fo{k) (43) 


iVg iV^ 

(k) ^ ^ -a ^ Wg (44) 

1=1 1=1 


Thus the power spectrum multipoles estimators for £ = 
0, 2,4 are 


h{k) = ^ 

4 22 


dVLk 


|Fo(k)p-iVo 


(45) 
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^2(fc) = I ^ F2{k)F*{k) (46) 

hik) = A I ^ F4(k)F*(k) (47) 

Ab(fc) = ^/ ^ |i^2(k)p-P2(fc)-^Po(fc) (48) 
where the standard normalization constant is m 



Nr. 


I 22 EE a^n(xj)^ 
i=i 


(49) 


and the shot noise obtained from the self-pairs in the first 
term in Eq. (45) 


Ng 

i=i i=i 

with the first term representing the true shot noise of 
galaxies, the second that of random objects. Replac¬ 
ing the true galaxy shot noise by its expectation value 

*0 ^0 = a(l + a)S^i the 

standard result m- However, using the true shot noise 
should always be preferred as it uses more information 
about the data [22]. Note that in general one would have 
a shot noise 




(51) 


i=i 


1=1 


FIG. 1: The bias ratio 64 = {Pib)/{P 4 ) between the two es¬ 
timators of the hexadecapole (symbols with error bars), and 
the ratio of their cosmic variance crl = {AP 41 ,)/{APi) as a 
function of k for the Las Damas LRG (Mg < —21.8) DR7 
mock catalogs. 


2.0 


5: " 


o; ' 

5: 1 



°'8.00 0.05 0.10 0.15 0.20 0.25 0.30 


k [h/Mpc] 


but this vanishes for £ > 0 . 

The implementation of the above sums by FFTs is 
straightforward and follows standard practice, e.g. the 
objects (galaxies or random) are interpolated to a grid 
to obtain a number density estimator n at gridpoints x, 
so that for any (tensor) quantity T 


Ng 

(E-E) 


1=1 1=1 


Tj 




— n(x) T(x) e' 


ik-: 


(52) 

where n = aur and the sum over the gridpoints x is 
performed using FFTs. In our implementation we use 
fourth-order interpolation to interlaced grids, which has 
superb antialiasing properties [53] . Compared to simpler 
interpolations, e.g. second-order (cloud in cell), fourth- 
order interpolation costs a factor of 8 more in computa¬ 
tional time and interlacing another factor of 2 , but this 
allows us to go up to the Nyquist frequency without any 
significant bias (thus saving a factor of at least eight due 
to the reduced size of the FFT to reach the same physical 
k). A detailed analysis of interpolation techniques and 


FIG. 2: Sames as Fig.l^for the PTHalos GMASS DRll mock 
catalogs. 


their impact on clustering properties will be presented 
elsewhere [24] . 

We now compare the two hexadecapole estimators P 4 
and Pih in terms of their expectation value and cosmic 
variance. For this purpose we use two sets of mock cata¬ 
logs: the Las Damas LRG {Mg < —21.8) DR7 mocks 
catalogs^ (160 realizations) and the PTHalos CMASS 
DRll mocks catalogs^ (600 realizations). In both cases 
we only use their north galactic cup versions. Figures [^ 
and [ 2 ] shows the results. We see that while both estima¬ 
tors agree within the errors on their expectation value, 


^ publicly available at http://lss.phy.vanderbilt.edu/lasdamas 
^ publicly available at http://www.marcmanera.net/mocks. We use 
version 11.1, see m for more details. 




































7 


the cosmic variance of P 4 is smaller than for Pjf, with the 
difference being smaller for the higher redshift and more 
dense sample (CMASS). Therefore the more computa¬ 
tionally expensive P 4 is preferred over the cheaper P 4 f,. 
This is not a very significant shortcoming as P 4 , is still 
orders of magnitude faster than traditional iV^ estimates. 


B. Bispectrum 


The bispectrum multipoles estimator is given by, 
3 

5(0) 


Bfi =n 


%i^Po(qi)Po(q2)Po(q3) - 

(53) 


-^ 1^3 -^33 


where the shot noise term is given by, 
r(0) _ TT /" ^3„ <5D(qi23) " 


= n / 


.•_i ^ ki 


NI 23 h3 L 


-Fb(qi)Po“'(-qi) + cyc. 


2 

ha 


Ng AT, 




i=i 

and for higher multipoles we obtain, 
3 


B[g = (2£+l)[] / 


3 ^ <5D(qi23) 


<23 ha 


-fVi 


w 

123 


(54) 


^€(qi)Po(q2)Po(q3) 
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with the shot noise given by (£ > 0 ), 
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where 

<(91, q) = 


Na 


Nr 


w'j Ci{qi ■ Xj) (57) 




i=i 


and ^(gi,qi) = P/'(qi). If desired the estimator in 
Eq. (55) can be symmetrized over its arguments in the 
obvious way. In the plane-parallel limit, this estimator 
for .£ = 2 reduces to the one in [ 20 ], which was used 
to measure the bispectrum quadrupole in Nbody simula¬ 
tions. 

We can simplify Eq. (56) by assuming the thin-shell 
approximation, which leads to 

d^qi P^(qi)Po“(-qi) 


Nig (2^+1) /■ ^ 
Jki Afej 
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FIG. 3: Reduced bispectrum multipoles Q^ig {i = 0, 2,4 from 
top to bottom) for galaxies in the LasDamas LRG (Mg < 
— 21.8) DR7 mock catalogs. The triangles correspond to ki = 
0.047 hMpc”^, A :2 = 2 fci as a function of the angle 9 between 
ki and k2. 


In computing the bispectrum, an additional complex¬ 
ity over the power spectrum case is that a search over 
closed triangles has to be done m, enforced by the 
delta function in Eqs. (53) and (55). There are many 


ways to do this, let us briefly discuss two options that 
we have found reasonably efhcient and implemented in 
the past |?7h29j . One is to use a fast (iVlniV) algorithm 
such as quicksort to sort Fourier coefficients into shells to 
quickly find qa = —qi 2 in shell h given qi in shell ki 
and q 2 in shell fc 2 . The sorting can be done once and the 
results stored in disk when multiple realizations need to 
be run (e.g. when running on mock catalogs) since it only 
depends on grid and bin size. The other is to use FFTs 
themselves to find closed triangles by using the plane- 
wave representation of the delta function and factorizing 
the estimator in real space, i.e. 


/ d^9*<5D(qi23)^>i(qi)^>2(q2)^4(q3) 

j_l Jki 

= <'’(x)F<«(x)f-l«(x) (59) 


where 


Fg\K)= [ d^qF^-^Fe{(l) (60) 

Jk 

and thus for each bin ki one must do an inverse FFT to 
(t) 

find the F^ , and then sum over real space to obtain the 
bispectrum for a given ki,k 2 ,k 3 . This estimator can be 





















trivially extended to hie:her-order spectra, as the trispec¬ 
trum [30]. 

Figure shows the results of measuring the bispec¬ 
trum for .^ = 0,2,4 for the LasDamas LRG {Mg < —21.8) 
DR7 mock catalogs for triangles that correspond to ki = 
0.047 h Mpc“^, fc 2 = 2 fci as a function of the angle 9 be¬ 
tween ki and k 2 . The bispectrum multipoles were com¬ 
puted using Eq. ( |^ with the shot noise correction given 
by Eq. (58). We see that the typical dependence on tri¬ 
angle shape of the reduced bispectrum defined as 


pC^) 

nW __ '°i23 _ 

Po(fci)Po(fc2) + Po{k2)Poik3) + Po{k3)Po{ki) 

( 61 ) 

is shared among the three multipoles. This was already 
known up to ^ = 2 in the plane-parallel approxima¬ 
tion . We leave the comparison of these measurements 
to theoretical predictions for an upcoming paper. 


C. Performance 


To get an idea of performance of the estimators pre¬ 
sented here, we now discuss for definiteness run times 
of our codes for the northern part of the CMASS sam¬ 
ple in the Baryon Oscillation Spectroscopic Survey ED, 
which comprises about Ng = 650,000 galaxies (and 
Nr « 100A"g). All timings in what follows are for a 
single processing core. The galaxies and random ob¬ 
jects are put into a cartesian box of 3.6 Gpc a side 
and fourth-order interpolated into interlaced 360^ grids, 
enough to reach k = 0.3 h Mpc~^ for the power spectrum 
and k = 0.2 /iMpc”^ for the bispectrum. To give an idea 
of how far we are from the plane-parallel approximation, 
for the DRll mask we have 


^xx — v.uvy, Lyy ~ 0.275, Lzz — 0.222, 
L,y~0.02, Lj,,--0.025, L,, ~-0.24, (62) 


where 


'' n '•h 

3 


ik-Xj 


Lao ^ 


‘'j 


' Wn e 


ik-: 


(63) 


which shows dominance of the x-direction but significant 
amplitudes in the other directions as well. 

The galaxies and randoms are interpolated and the 7 
EFTs computed. For the galaxies this takes 22 seconds, 
while for the randoms 210 seconds. Computing the addi¬ 
tional 15 FFTs (and interpolations) as well takes about 
60 seconds in total for the galaxies, while computing the 
corresponding FFTs for ^“(q) (needed for bispectrum 
shot noise subtraction) doubles the timings again (for 
a total of 44 interpolations and FFTs). The overden¬ 
sity fields are constructed and the £ = 0,2,4 multipoles 
from fcniin = 0.0052 hMpc”^ up to fcmax = 0.3/iMpc“^ 
in Sk = 0.0052 h Mpc“^ bins computed in 2 seconds. For 
the bispectrum we compute monopole and quadrupole for 


all triangles with sides between fcmin = 0.0052 hMpc“^ 
and /cmax = 0.209 h Mpc“^ in 5k = 0.0052 h Mpc“^ bins 
in about 11 minutes; adding the hexadecapole yields an 
additional 5 minutes. Clearly, this is a more than ad¬ 
equate performance and can be scaled to significantly 
larger data sets with no foreseeable issues (in fact we 
have routinely used these algorithms for tens of billion 
particle simulations in the plane-parallel case). In ad¬ 
dition, most importantly, this performance allows us to 
estimate the covariance matrix of these estimators using 
~ 10^ survey realizations with realistic masks and noise 
properties, as has been already presented in the past for 
smaller surveys EZlEa. 


IV. CONCLUSIONS 

Computing the redshift-space clustering of galaxies is 
the key goal of large-scale structure and one of the most 
sensitive probes to test gravity and dark energy, bias and 
primordial non-Gaussianity. Application to wide surveys 
demands defining statistics that characterize the effect of 
redshift-space distortions allowing for the angular depen¬ 
dence of the line of sight along which distortions operate. 

In this paper, starting from a definition of local estima¬ 
tors that generalize the dehnition of power spectrum and 
bispectrum to the case of lack of statistical homogeneity, 
we defined natural estimators for power spectrum and 
bispectrum multipoles. In the case of the power spec¬ 
trum multipoles, our natural estimator agrees for i = 2 
with [16] . while for the bispectrum multipoles it is new. 
We then considered slight modifications to these estima¬ 
tors in which the “center of mass line of sight” is allowed 
to rotate among the different members of a pair (for the 
power spectrum) or triplet (for the bispectrum). As a 
result of this, we presented very efficient estimators to 
calculate the power spectrum and bispectrum multipoles, 
which require only FFTs. 

For the power spectrum quadrupole, an additional 6 
FFTs are required over the monopole, while for the power 
spectrum hexadecapole we presented two estimators, one 
that requires an additional 15 FFTs over the quadrupole, 
and another one which does not require any additional 
FFTs. We implemented both in mock catalogs with 
realistic geometries and found that while the two hex¬ 
adecapole estimators agree with each other in the mean 
value, the more expensive estimator has a somewhat 
lower cosmic variance, and thus higher signal to noise. 
We also showed first results for the bispectrum £ = 0, 2,4 
estimators in mock catalogs and presented specihcs about 
the performance of such estimators for the BOSS sur¬ 
vey. Their speed and scaling makes application to larger 
future datasets such as eBOSS, Euclid and DESI quite 
encouraging. In the near future we will present detailed 
studies of the bias of these estimators when compared to 
the plane-parallel approximation (which is always used to 
make theoretical predictions) as well as their application 
to the DR12 sample of BOSS galaxies. 
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When this paper was in preparation |33] appeared on 
the arXiv with the same results for the FFT factorization 
of the power spectrum multipoles estimator leading to 
our Eqs. (201 and (23). 
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